import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import random

# Generate a random number between a and b
def GetRandomFactor():
    random.seed()  # Initialize the random number generator using system time
    return random.uniform(0.9990, 1.0001)

def GetRandomNumber():
    random.seed()  # Initialize the random number generator using system time
    return random.uniform(-0.1, 0.1)

plt.rcParams['font.sans-serif'] = ['SimHei']  # Use SimHei font
plt.rcParams['font.weight'] = 'bold'  # Set font weight to bold
plt.rcParams['axes.unicode_minus'] = False  # Solve the problem of displaying minus sign

# Define the season function
def season_factor(t):
    # Divide the year into four seasons, each lasting 3 months
    season = (t % 12) // 3  # Each season lasts 3 months
    return season

# Define parameters that change with time
def r(t):
    season = season_factor(t)

    if season == 1:
        increase_rate = 0.2
    elif season == 2:
        increase_rate = 0.24
    elif season == 3:
        increase_rate = 0.25
    else:
        increase_rate = 0.16
    return 0.3 + increase_rate * GetRandomFactor() * 0.7 + 1.2 * np.sin(0.005 * t)

# Define the system of differential equations
def evolution_forests_to_farmland(y0, t, params):
    P, I, B, S = y0  # Set initial population sizes
    K, a, C, a_P, b, I_M, yp_I, rr, yp, Ba, B_M, mu_B, k, S_M, lam = params

    '''
    K: Carrying capacity of plants
    a: Pesticide impact factor on plants
    C: Chemical concentration
    a_P: Plant resistance to insects
    b: Insect growth rate
    I_M: Maximum carrying capacity of insects
    yp_I: Insect natural mortality rate
    rr: Chemical impact factor on insects
    yp: Bird predation efficiency on insects
    Ba: Bird growth rate
    B_M: Maximum carrying capacity of birds
    mu_B: Bird natural mortality rate
    k: Soil recovery rate
    S_M: Maximum soil fertility
    lam: Chemical impact factor on soil
    '''
    dPdt = r(t) * P * (1 - P / K) - a * P * C - a_P * P * I  # Plant growth
    dIdt = b * I * (1 - I / I_M) - rr * I * C - yp * B * I - yp_I * I + (P / 300) * 15  # Insect growth
    dBdt = Ba * B * (1 - B / B_M) + yp * B * I - mu_B * B  # Bird growth
    dSdt = k * (S_M - S) - lam * S * C  # Soil quality

    return [dPdt, dIdt, dBdt, dSdt]

# Initial population sizes
P0 = 350  # Plants
I0 = 130  # Insects
B0 = 20  # Birds
S0 = 1.2  # Soil fertility
y0 = [P0, I0, B0, S0]

# Time points
t = np.linspace(0, 120, 1000)

# Parameter settings
params = (
    500,  # K: Carrying capacity of plants
    3,  # a: Pesticide impact factor on plants
    0.01,  # C: Chemical concentration
    0.0008,  # a_P: Plant resistance to insects
    0.3,  # b: Insect growth rate
    200,  # I_M: Maximum carrying capacity of insects
    0.05,  # yp_I: Insect natural mortality rate
    20.8,  # rr: Chemical impact factor on insects
    0.00005,  # yp: Bird predation efficiency on insects(too many)
    0.1,  # Ba: Bird growth rate
    50,  # B_M: Maximum carrying capacity of birds
    0.06,  # mu_B: Bird natural mortality rate
    0.13,  # k: Soil recovery rate
    1.5,  # S_M: Maximum soil fertility
    2.8  # lam: Chemical impact factor on soil
)

# Solve the system of differential equations
solution = odeint(evolution_forests_to_farmland, y0, t, args=(params,))

# Plot the results
plt.figure(figsize=(10, 8))
plt.plot(t, solution[:, 0] / P0, label='Plant current biomass / initial biomass (P)', color='blue')
plt.plot(t, (solution[:, 1] / I0) * GetRandomFactor() + GetRandomNumber() * 0.55 * GetRandomFactor() * np.sin(t), label='Insect current biomass / initial biomass (I)', color='orange')
plt.plot(t, solution[:, 2] / B0 + GetRandomNumber() * 0.05 * GetRandomFactor() * np.sin(t), label='Bird(bats) current biomass / initial biomass (B)', color='green')
plt.plot(t, solution[:, 3] / S0, label='Soil fertility / initial soil fertility (S)', color='red')

plt.legend()

plt.xlabel('Time (months)', fontweight='bold')
plt.ylabel('Initial biomass (soil fertility) / current biomass (soil fertility)', fontweight='bold')
plt.title('Population size and soil fertility over time in a mature agroecosystem (using chemicals)', fontweight='bold')

# Set y-axis ticks to 0.1 intervals
plt.yticks(np.arange(0.8, 1.5, 0.1))

# Set x-axis ticks to 12 intervals
plt.xticks(np.arange(0, 121, 12))

plt.show()